Twin analysis using TwinAnalysis

Ivan Voronin

12 января, 2021

Contents

Twin method: Fundamentals

ACE

ADE

Univariate ACE

Univariate ACE

Univariate ACE

Univariate ACE

Univariate ACE

TwinAnalysis

Univariate toy data

# Making univariate toy data
library(dplyr,
        quietly = TRUE)

mz_data <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0),
  Sigma = matrix(c(1.0, 0.8,
                   0.8, 1.0),
                 nrow = 2)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'X2'))

dz_data <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0),
  Sigma = matrix(c(1.0, 0.5,
                   0.5, 1.0),
                 nrow = 2)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'X2'))

data <- rbind(data.frame(zyg = 'MZ', mz_data),
              data.frame(zyg = 'DZ', dz_data)) %>%
  mutate(
    zyg = factor(zyg, levels = c('MZ', 'DZ'))
  )

head(data)
##   zyg          X1         X2
## 1  MZ  1.97640719  0.8614917
## 2  MZ  0.19813515  0.5722313
## 3  MZ -0.21476752 -0.1499316
## 4  MZ -0.05976396 -1.6416417
## 5  MZ -0.10528725 -0.8397834
## 6  MZ  1.39832330  1.1882021

Descriptives

# Installing the packages
# install.packages('devtools')
# library(devtools)
# install_github('ivanvoronin/mlth.data.frame')
# install_github('ivanvoronin/TwinAnalysis')

# Loading the packages
library(OpenMx)
library(TwinAnalysis)

# Descriptive statistics
get_univ_descriptives(data, zyg = 'zyg', vars = 'X')
## $X
##    Twin1------------------ Twin2------------------- Cor--------------------------
##    n   M          SD       n   M          SD        r         lowerCI   upperCI  
## MZ 100  0.2892160 0.996696 100  0.1193898 0.8943036 0.7601622 0.6626901 0.8323086
## DZ 100 -0.1712873 1.067547 100 -0.1975880 0.9445352 0.4553250 0.2843362 0.5982401

Univariate ACE

# Univariate ACE model
model <- univariate_ace(data, zyg = 'zyg', ph = 'X')

model <- mxRun(model)

summary(model)
## Summary of ACE 
##  
## free parameters:
##            name matrix row col     Estimate  Std.Error A
## 1 ACE.mean[1,1]   mean   1   1 -0.008753464 0.06321116  
## 2    ACE.a[1,1]      a   X   X  0.785090384 0.10184093  
## 3    ACE.c[1,1]      c   X   X  0.381478731 0.19748414  
## 4    ACE.e[1,1]      e   X   X  0.482854666 0.03427370  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:              4                    396
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:              1021.598
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/400
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       229.5979               1029.598
## BIC:     -1076.5358               1042.791
##       |  Sample-Size Adjusted
## AIC:                 1029.803
## BIC:                 1030.119
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:43:25 
## Wall clock time: 0.07648754 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Output - parameters

# Output tables from the univariate ACE
get_output_tables(model)
## $Variance_components
##        A(%)      C(%)      E(%)
## X 0.6194384 0.1462512 0.2343104
## 
## $Raw_variance
##       Total         A        C         E
## X 0.9950416 0.6163669 0.145526 0.2331486

Output - model fit

# Model fit
model <- ref_models(model)

fit_stats(model)
##   model      base ep minus2LL  df      AIC       BIC
## 1   ACE Saturated  4 1021.598 396 229.5979 -1076.536
##        CFI      TLI     RMSEA   diffLL diffdf           p
## 1 0.857794 0.952598 0.1128539 21.28321      6 0.001631543

Constrained models

# Constrained models
# This is a list of MxModels
nested_models <- def_nested_models(
  model,
  AE = mxConstraint(c[1, 1] == 0),
  CE = mxConstraint(a[1, 1] == 0),
  E = list(mxConstraint(c[1, 1] == 0),
           mxConstraint(a[1, 1] == 0)),
  run = TRUE
)

fit_stats(model, nested_models = nested_models)
## Warning in computeFitStatistics(likelihood, DoF, chi,
## chiDoF, retval[["numObs"]], : Your model may be mis-
## specified (and fit worse than an independence model),
## or you may be using the wrong independence model, see ?
## mxRefModels
##   model      base ep minus2LL  df      AIC        BIC
## 1   ACE Saturated  4 1021.598 396 229.5979 -1076.5358
## 2    AE       ACE  4 1022.470 397 228.4697 -1080.9623
## 3    CE       ACE  4 1039.653 397 245.6534 -1063.7786
## 4     E       ACE  4 1130.385 398 334.3846  -978.3457
##          CFI       TLI     RMSEA      diffLL diffdf
## 1  0.8577940 0.9525980 0.1128539  21.2832071      6
## 2  0.8589871 0.9597106 0.1040432   0.8717727      1
## 3  0.6990974 0.9140278 0.1519838  18.0554872      1
## 4 -0.1358265 0.7160434 0.2762131 108.7866932      2
##              p
## 1 1.631543e-03
## 2 3.504650e-01
## 3 2.145594e-05
## 4 2.383799e-24

Bivariate toy data

# Making bivariate toy data
mz_data2 <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0, 0, 0),
  Sigma = matrix(
    c(1.0, 0.7, 0.6, 0.5,
      0.7, 1.0, 0.5, 0.8,
      0.6, 0.5, 1.0, 0.7,
      0.5, 0.8, 0.7, 1.0),
    nrow = 4)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'Y1', 'X2', 'Y2'))

dz_data2 <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0, 0, 0),
  Sigma = matrix(
    c(1.0, 0.7, 0.4, 0.2,
      0.7, 1.0, 0.2, 0.4,
      0.4, 0.2, 1.0, 0.7,
      0.2, 0.4, 0.7, 1.0),
    nrow = 4)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'Y1', 'X2', 'Y2'))

data2 <- rbind(data.frame(zyg = 'MZ', mz_data2),
               data.frame(zyg = 'DZ', dz_data2)) %>%
  mutate(
    zyg = factor(zyg, levels = c('MZ', 'DZ'))
  )

head(data2)
##   zyg         X1          Y1          X2          Y2
## 1  MZ -0.5435463 -0.41234556 -0.86241804 -0.29155328
## 2  MZ  0.3948285 -0.84486187 -0.18510648 -1.11259359
## 3  MZ  0.3981069 -0.08203682 -0.05917079 -0.18544979
## 4  MZ -0.6577894 -0.73082395 -0.02673546  0.03067235
## 5  MZ  0.4123762  0.71153558  0.41007149  0.53892970
## 6  MZ -0.1671188 -1.43100131 -0.50763096 -1.31595716

Descriptives

# Descriptive statistics
get_univ_descriptives(data2, 
                      zyg = 'zyg', 
                      vars = c('X', 'Y'))
## $X
##    Twin1------------------- Twin2-------------------- Cor--------------------------
##    n   M          SD        n   M           SD        r         lowerCI   upperCI  
## MZ 100  0.1602093 0.9891965 100  0.09714586 0.9135136 0.5967201 0.4534493 0.7099298
## DZ 100 -0.1450704 0.9265116 100 -0.20180809 1.0229349 0.4750329 0.3072861 0.6141479
## 
## $Y
##    Twin1------------------ Twin2-------------------- Cor--------------------------
##    n   M          SD       n   M           SD        r         lowerCI   upperCI  
## MZ 100  0.1208040 1.101084 100  0.06408469 1.0693377 0.8602316 0.7987792 0.9039190
## DZ 100 -0.1001922 0.983983 100 -0.24580522 0.9595227 0.4185444 0.2420229 0.5682471

Twin correlations

# Cross-twin cross-trait correlations
get_cross_trait_cors(data2,
                     zyg = 'zyg',
                     vars = c('X', 'Y'))
## zyg: MZ
##           X1        Y1        X2        Y2
## X1 1.0000000 0.6625874 0.5967201 0.4859931
## Y1 0.6625874 1.0000000 0.5098315 0.8602316
## X2 0.5967201 0.5098315 1.0000000 0.6378926
## Y2 0.4859931 0.8602316 0.6378926 1.0000000
## -------------------------------------------- 
## zyg: DZ
##           X1        Y1        X2        Y2
## X1 1.0000000 0.6686251 0.4750329 0.2190231
## Y1 0.6686251 1.0000000 0.2353515 0.4185444
## X2 0.4750329 0.2353515 1.0000000 0.6493918
## Y2 0.2190231 0.4185444 0.6493918 1.0000000

Bivariate ACE

# Bivariate ACE model
model2 <- multivariate_ace(
  data = data2,
  zyg = 'zyg',
  ph = c('X', 'Y')
)

model2 <- mxRun(model2)

summary(model2)
## Summary of ACE 
##  
## free parameters:
##             name matrix row col      Estimate  Std.Error A
## 1  ACE.mean[1,1]   mean   1   1 -3.293302e-02 0.06011835  
## 2  ACE.mean[1,2]   mean   1   2 -5.684037e-02 0.06549453  
## 3     ACE.a[1,1]      a   X   X  2.041926e-07 0.24116291  
## 4     ACE.a[2,1]      a   Y   X  5.774286e-01 0.10958100  
## 5     ACE.a[2,2]      a   Y   Y  9.063612e-01 0.09302743  
## 6     ACE.c[1,1]      c   X   X  4.986504e-01 0.24273142  
## 7     ACE.c[2,1]      c   Y   X -9.775818e-02 0.63518951  
## 8     ACE.c[2,2]      c   Y   Y  2.565577e-01 0.32877565  
## 9     ACE.e[1,1]      e   X   X  4.475017e-01 0.02362679  
## 10    ACE.e[2,1]      e   Y   X  3.870517e-01 0.04866025  
## 11    ACE.e[2,2]      e   Y   Y  4.042446e-01 0.02815001  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:             11                    789
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:               1793.14
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/800
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       215.1395               1815.140
## BIC:     -2387.2329               1851.421
##       |  Sample-Size Adjusted
## AIC:                 1816.544
## BIC:                 1816.572
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:43:25 
## Wall clock time: 0.07370257 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Output - parameters

# Output tables from the bivariate ACE
get_output_tables(model2)
## $Variance_components
##        A(%)       C(%)      E(%)
## X 0.3540660 0.27419454 0.3717394
## Y 0.7818313 0.06264417 0.1555245
## 
## $Covariation_components
##             Total         A           C         E
## X <-> X 0.9416995 0.3334238  0.25820885 0.3500668
## X <-> Y 0.6547418 0.5233589 -0.02508062 0.1564636
## Y <-> Y 1.0507262 0.8214906  0.06582187 0.1634137
##              A(%)       C(%)      E(%)
## X <-> X 0.3540660 0.27419454 0.3717394
## X <-> Y 0.7424551 0.03558023 0.2219647
## Y <-> Y 0.7818313 0.06264417 0.1555245
## 
## $All_correlations
##         r_A        r_C       r_E     Total
## X <-> X   1  1.0000000 1.0000000 1.0000000
## X <-> Y   1 -0.1923834 0.6541744 0.6582171
## Y <-> Y   1  1.0000000 1.0000000 1.0000000

Output - model fit

# Model fit
model2 <- ref_models(model2)

fit_stats(model2)
##   model      base ep minus2LL  df      AIC       BIC
## 1   ACE Saturated 11  1793.14 789 215.1395 -2387.233
##         CFI       TLI      RMSEA   diffLL diffdf       p
## 1 0.9970143 0.9978924 0.02022251 18.39043     17 0.36462

Vignette

OpenMx

A model (mxModel):

  • matrices (mxMatrix)
  • matrix algebra (mxAlgebra)
  • data (mxData)
  • fit function (mxFitFunction…)
  • expectation (mxExpectation…)
  • sub-models (mxModel)
  • CIs (mxCI), constraints (mxConstraint)

Univatiate ACE

# Twin models in OpenMx
library(OpenMx)

# Univariate ACE model
model <- mxModel(
  name = 'ACE',
  
  # Means
  mxMatrix(type = 'Full',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'mean'),
  mxAlgebra(cbind(mean, mean),
            name = 'expMeans'),
  
  # Variance components
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'a'),
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'c'),
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'e'),
  mxAlgebra(t(a) %*% a, name = 'A'),
  mxAlgebra(t(c) %*% c, name = 'C'),
  mxAlgebra(t(e) %*% e, name = 'E'),
  
  # Expected covariance
  mxAlgebra(rbind(cbind(A + C + E, A + C    ),
                  cbind(A + C,     A + C + E)),
            name = 'expCovMZ'),
  mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
                  cbind(.5%x%A + C, A + C + E)),
            name = 'expCovDZ'),
  
  # MZ submodel
  mxModel(name = 'MZ',
          mxData(observed = mz_data,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovMZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'X2')),
          mxFitFunctionML()),
  
  # DZ submodel
  mxModel(name = 'DZ',
          mxData(observed = dz_data,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovDZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'X2')),
          mxFitFunctionML()),
  
  mxFitFunctionMultigroup(c('MZ','DZ')),
  
  # Output tables
  mxAlgebra(A + C + E, name = 'V'),
  mxAlgebra(cbind(V, A, C, E),
            dimnames = list('X', c('Total', 'A', 'C', 'E')),
            name = 'Raw_variance'),
  mxAlgebra(cbind(A / V, C / V, E / V),
            dimnames = list('X', c('A(%)', 'C(%)', 'E(%)')),
            name = 'Variance_components')
)

model <- mxRun(model)

summary(model)
## Summary of ACE 
##  
## free parameters:
##            name matrix row col     Estimate  Std.Error A
## 1 ACE.mean[1,1]   mean   1   1 -0.008753687 0.06321119  
## 2    ACE.a[1,1]      a   1   1  0.785090556 0.10183866  
## 3    ACE.c[1,1]      c   1   1  0.381478510 0.19747930  
## 4    ACE.e[1,1]      e   1   1  0.482854720 0.03427365  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:              4                    396
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:              1021.598
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/400
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       229.5979               1029.598
## BIC:     -1076.5358               1042.791
##       |  Sample-Size Adjusted
## AIC:                 1029.803
## BIC:                 1030.119
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:43:25 
## Wall clock time: 0.04036498 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Univariate output

# Univariate output
mxEval(Raw_variance, model)
##       Total         A         C         E
## X 0.9950417 0.6163672 0.1455259 0.2331487
mxEval(Variance_components, model)
##        A(%)     C(%)      E(%)
## X 0.6194385 0.146251 0.2343105

Bivariate ACE

# Bivariate ACE model
model2 <- mxModel(
  name = 'ACE',
  
  # Means
  mxMatrix(type = 'Full',
           nrow = 1, ncol = 2,
           free = TRUE,
           name = 'mean'),
  mxAlgebra(cbind(mean, mean),
            name = 'expMeans'),
  
  # Variance components
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'a'),
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'c'),
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'e'),
  mxAlgebra(t(a) %*% a, name = 'A'),
  mxAlgebra(t(c) %*% c, name = 'C'),
  mxAlgebra(t(e) %*% e, name = 'E'),
  
  # Covariance
  mxAlgebra(rbind(cbind(A + C + E, A + C    ),
                  cbind(A + C,     A + C + E)),
            name = 'expCovMZ'),
  mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
                  cbind(.5%x%A + C, A + C + E)),
            name = 'expCovDZ'),
  
  # MZ submodel
  mxModel(name = 'MZ',
          mxData(observed = mz_data2,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovMZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'Y1', 'X2', 'Y2')),
          mxFitFunctionML()),
  
  # DZ submodel
  mxModel(name = 'DZ',
          mxData(observed = dz_data2,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovDZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'Y1', 'X2', 'Y2')),
          mxFitFunctionML()),
  
  mxFitFunctionMultigroup(c('MZ','DZ')),
  
  # Output tables
  mxAlgebra(A + C + E, name = 'V'),
  mxAlgebra(abs(A) + abs(C) + abs(E), name = 'Vabs'),
  mxAlgebra(abs(A) / Vabs, name = 'Aperc'),
  mxAlgebra(abs(C) / Vabs, name = 'Cperc'),
  mxAlgebra(abs(E) / Vabs, name = 'Eperc'),
  mxAlgebra(cbind(diag2vec(Aperc),
                  diag2vec(Cperc),
                  diag2vec(Eperc)),
            dimnames = list(c('X', 'Y'), c('A(%)', 'C(%)', 'E(%)')),
            name = 'Variance_components'),
  
  mxAlgebra(cbind(vech(V),
                  vech(A), vech(C), vech(E),
                  vech(Aperc), vech(Cperc), vech(Eperc)),
            dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
                            c('Total', 'A', 'C', 'E', 'A(%)', 'C(%)', 'E(%)')),
            name = 'Covariation_components'),
  
  mxMatrix(type = "Iden",
           nrow = 2, ncol = 2,
           name = "I"),
  mxAlgebra(sqrt(I * V), name = "SD_total"),
  mxAlgebra(sqrt(I * A), name = 'SD_A'),
  mxAlgebra(sqrt(I * C), name = 'SD_C'),
  mxAlgebra(sqrt(I * E), name = 'SD_E'),
  
  mxAlgebra(solve(SD_total) %&% V, name='r_total'),
  mxAlgebra(solve(SD_A) %&% A, name = 'r_A'),
  mxAlgebra(solve(SD_C) %&% C, name = 'r_C'),
  mxAlgebra(solve(SD_E) %&% E, name = 'r_E'),
  
  mxAlgebra(cbind(vech(r_A), vech(r_C), vech(r_E),
                  vech(r_total)),
            dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
                            c('r_A', 'r_C', 'r_E', 'Total')),
            name='All_correlations')
)

model2 <- mxRun(model2)

summary(model2)
## Summary of ACE 
##  
## free parameters:
##             name matrix row col      Estimate  Std.Error A
## 1  ACE.mean[1,1]   mean   1   1 -3.293282e-02 0.06011839  
## 2  ACE.mean[1,2]   mean   1   2 -5.684027e-02 0.06549455  
## 3     ACE.a[1,1]      a   1   1 -1.357998e-07 0.24112224  
## 4     ACE.a[2,1]      a   2   1  5.774276e-01 0.10955114  
## 5     ACE.a[2,2]      a   2   2  9.063608e-01 0.09299972  
## 6     ACE.c[1,1]      c   1   1  4.986521e-01 0.24264231  
## 7     ACE.c[2,1]      c   2   1 -9.775422e-02 0.63495296  
## 8     ACE.c[2,2]      c   2   2  2.565589e-01 0.32867213  
## 9     ACE.e[1,1]      e   1   1  4.475017e-01 0.02362661  
## 10    ACE.e[2,1]      e   2   1  3.870517e-01 0.04865963  
## 11    ACE.e[2,2]      e   2   2  4.042445e-01 0.02814985  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:             11                    789
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:               1793.14
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/800
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       215.1395               1815.140
## BIC:     -2387.2329               1851.421
##       |  Sample-Size Adjusted
## AIC:                 1816.544
## BIC:                 1816.572
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:43:26 
## Wall clock time: 0.06721616 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Bivariate output

# Bivariate output
mxEval(Variance_components, model2)
##        A(%)       C(%)      E(%)
## X 0.3540649 0.27419561 0.3717395
## Y 0.7818308 0.06264473 0.1555245
mxEval(All_correlations, model2)
##         r_A        r_C       r_E    Total
## X <-> X   1  1.0000000 1.0000000 1.000000
## X <-> Y   1 -0.1923752 0.6541743 0.658217
## Y <-> Y   1  1.0000000 1.0000000 1.000000

Bivariate output

mxEval(Covariation_components, model2)
##             Total         A           C         E
## X <-> X 0.9416993 0.3334227  0.25820980 0.3500668
## X <-> Y 0.6547416 0.5233578 -0.02507971 0.1564635
## Y <-> Y 1.0507260 0.8214899  0.06582245 0.1634136
##              A(%)       C(%)      E(%)
## X <-> X 0.3540649 0.27419561 0.3717395
## X <-> Y 0.7424557 0.03557906 0.2219653
## Y <-> Y 0.7818308 0.06264473 0.1555245

mlth.data.frame

library(mlth.data.frame)

# mlth.data.frame - multiheaded data.frame
model2 <- multivariate_ace(
  data = data2,
  zyg = 'zyg',
  ph = c('X', 'Y')
)
model2 <- def_ci(model2, ci = 'Variance_components')
model2 <- mxRun(model2, intervals = TRUE)

get_output_tables(model2)
## $Variance_components
##   A(%)------------------------- C(%)----------------------------- E(%)-------------------------
##   est       lCI       uCI       est        lCI          uCI       est       lCI       uCI      
## X 0.3540660 0.1388926 0.5823833 0.27419454 9.228158e-02 0.4861962 0.3717394 0.2892442 0.4715775
## Y 0.7818313 0.5048375 0.8786609 0.06264417 6.563101e-31 0.3347188 0.1555245 0.1145364 0.2128753
## 
## $Covariation_components
##             Total         A           C         E
## X <-> X 0.9416995 0.3334238  0.25820885 0.3500668
## X <-> Y 0.6547418 0.5233589 -0.02508062 0.1564636
## Y <-> Y 1.0507262 0.8214906  0.06582187 0.1634137
##              A(%)       C(%)      E(%)
## X <-> X 0.3540660 0.27419454 0.3717394
## X <-> Y 0.7424551 0.03558023 0.2219647
## Y <-> Y 0.7818313 0.06264417 0.1555245
## 
## $All_correlations
##         r_A        r_C       r_E     Total
## X <-> X   1  1.0000000 1.0000000 1.0000000
## X <-> Y   1 -0.1923834 0.6541744 0.6582171
## Y <-> Y   1  1.0000000 1.0000000 1.0000000

M <- get_output_tables(model2)$Variance_components

M[['A(%)']]
##   est       lCI       uCI      
## X 0.3540660 0.1388926 0.5823833
## Y 0.7818313 0.5048375 0.8786609

# Formatting the table for html or latex report
# Powered by knitr and kableExtra
library(kableExtra)

M %>% 
  behead() %>%
  kable(digits = 3,
        align = 'c') %>%
  add_complex_header_above(M)
A(%)
C(%)
E(%)
est lCI uCI est lCI uCI est lCI uCI
X 0.354 0.139 0.582 0.274 0.092 0.486 0.372 0.289 0.472
Y 0.782 0.505 0.879 0.063 0.000 0.335 0.156 0.115 0.213 

# Saving the table for the output
M %>%
  register_output(
    name = 'Table 1',
    caption = 'Table 1: Variance components')
##   A(%)------------------------- C(%)----------------------------- E(%)-------------------------
##   est       lCI       uCI       est        lCI          uCI       est       lCI       uCI      
## X 0.3540660 0.1388926 0.5823833 0.27419454 9.228158e-02 0.4861962 0.3717394 0.2892442 0.4715775
## Y 0.7818313 0.5048375 0.8786609 0.06264417 6.563101e-31 0.3347188 0.1555245 0.1145364 0.2128753

# Run at the end of the script
# to write all output tables into a single file
# openxlsx package required
write.xlsx.output(file = 'output.xlsx')

Vignettes